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Abstract 



Higher order nonlinear effects of site blocking, which induces oscillations in 
the Monte Carlo simulation of CO oxidation, are outlined here. It is shown 
that the rate equations which include these effects exhibit a supercritical 
Hopf bifurcation in parameter domains where Monte Carlo simulations lead 
to oscillations. 
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I. INTRODUCTION 



Surface chemical reactions play an important role in the processes of heterogeneous 
catalysis, which are widely used in the chemical industry [1]. Monte Carlo simulation is a 
powerful tool for a microscopic description of these reactions [2,3]. Kinetic phase transitions, 
oscillations and chaos, pattern formation, and coupling between catalytic oscillators have 
been studied by the Monte Carlo method [4-7]. 

The oxidation of carbon monoxide CO is one of the most extensively studied surface 
catalytic reactions [8,9]. In a recent Monte Carlo study, an attenpt of a mean field modeling 
was made for site blocking induced surface coverage oscillations [10]. The authors of Ref. [10] 
studied a mathematical model based on macroscopic rate equations for the characterization 
of oscillations seen in the Monte Carlo simulations. However, it was found that the mean 
field model does not lead to oscillatory dynamics in a parameter region where their Monte 
Carlo simulations show oscillations. In order to obtain oscillations, they numerically solved 
a part of their rate equations with data from the Monte Carlo simulation for dynamics of the 
free site and the inert molecule coverage. Such a combination of deterministic and stochastic 
time evolutions seems not to be satisfactory; it is more desirable to find a supercritical Hopf 
bifurcation in the model to ensure the occurrence and stability of oscillations. The aim of 
this work is to show that inclusion of high order nonlinear terms, which represent interactions 
inherent in this system, allows us to find oscillations within the rate equations. We find a 
Hopf bifurcation in the extended system, and we can determine the parameter domain for 
oscillatory dynamics. 



II. SITE BLOCKING INDUCED OSCILLATIONS 

The model system under study is the ZGB model [4], which is extended by the absorption 
and desorption of an inert molecule: 

co gas + * n CO ads , (1) 

CO ads CO gas + *, 

2 + 2* n 20 ads , 

CO ads + Oads ^ C0 2 + 2*, 
M gas + *^M ads , 
M ads ^ M gas + *, 

where ads indicates that the molecule is absorbed on the surface and gas indicates that the 
molecule is in the gas phase. * indicates a vacant site and M represents the site blocking 
inert molecule. 

As was reported by Jansen and Nieminen [10], one can observe sustained coverage oscil- 
lations in the Monte Carlo simulations of Eqn. (1). In Fig. 1 we show results of our Monte 
Carlo simulation for these coverage oscillations. All variables and time are dimensionless in 
this paper. As compared with Ref. [10], we have used a different set of parameters to show 
that these oscillations can be observed in wide range of parameters. To improve the numer- 
ical performance we have implemented lists in our Monte Carlo simulations. Separate lists 
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for the locations of free sites, CO, and inert molecules were updated during the simulations. 
For example, for an adsorption of a chosen adsorbent, a site was randomly selected from the 
list of free sites; for a desorption, a site was randomly selected from the corresponding list 
of CO or M. 

The authors of Ref. [10] studied rate equations for a macroscopic characterization of 
the oscillations seen in the Monte Carlo simulation. The rate equations for the coverage of 
CO ads - x, O ads - y, and M ads - z are [10]: 

a = l- x-y-z, (2) 
x = W\0 — w 2 x — Aw 3 xy, 
y = Aw^a 2 — Aw3xy, 
z = w 5 a — w 6 z, 

where a is the free site coverage. It is known that Eqn. (2) shows oscillations in a narrow 
parameter domain [11], however, temporal pattern of these oscillations are very different 
from the oscillations seen in Monte Carlo simulations as shown in Fig. 1. It was found that 
Eqn. (2) does not show oscillations in the parameter domain where Monte Carlo simulations 
lead to oscillations. To obtain oscillations, the authors of Ref. [10] used data for a and z 
from the Monte Carlo simulation. We suppose that the reason why Eqn. (2) does not 
show oscillations in the parameter domains of intereset is that it does not include important 
higher order nonlinear effects inherent in Eqn. (1), which are essential for oscillations shown 
in Fig. 1. 

III. HIGHER ORDER NONLINEAR EFFECTS 

Monte Carlo simulations of Eqn. (1) suggest that two nonlinear effects are crucial for 
oscillations in this system. They destabilize O and CO rich surfaces, correspondingly. 

In Fig. 2 we show snapshots of the surface at two different time moments. Fig. 2a 
corresponds to the time moment when there are no CO ads molecules left on the surface. 
Note that in the original ZGB model, for the parameters we have chosen, the system goes 
to the state fully covered by O ads [4]. However, the presence of the inert molecules shown 
in Fig. 2a as the brightest spots, prevent such a fully O ads covered state. Moreover, these 
inert molecules allow for the presence of enough vacant sites; these are shown as black areas 
in Fig. 2. In these vacant areas, CO molecules can be absorbed; their coverage can grow 
if they are isolated by the inert molecules from O ads covered areas. By taking into account 
this nonlinear effect, which includes interaction among CO ads , O ads , and M ads , the reaction 
term in Eqn. (2) can be modified as Awsxy(l — kz), where k is the blocking coefficient. 

Fig. 2b shows nucleation of the O ads covered state on the CO ads rich surface. As an 
adsorption of 2 requires two free sites, nucleation of the O ads island may seem an unlikely 
process, if one starts from the fully CO ads covered state as an initial condition. However, 
there is an another nonlinear effect which may lead to explosive growth of an O ads island. 
Suppose that the desorptions of two CO ads molecules have vacated a pair of adjacent sites. 
On these sites, 2 can be absorbed, and if this happens, CO ads and O ads will immediately 
react. As the surface is rich with CO ads , the reaction will free four sites. On these four sites, 
two 2 molecules can be absorbed. This cascade may lead to explosive growth of the O ads 
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island. Therefore, there is an another nonlinear effect which is particularly noticeable during 
a nucleation of an O a d s island on a CO a( i s rich surface. In the lowest order of nonlinearity, 
due to the nonlinear increase of available adjacent free sites through the reaction, the growth 
°f O a d s can be termed as 4u>w 3 x 2 y, where u is the coefficient for the nonlinear growth of 
adjacent free sites. We verified that Auw 3 x 2 y 2 can also serve well for this term. 

The two nonlinear effects manifest themselves randomly in Monte Carlo simulations. For 
a macrospcopic description, we consider them deterministic processes. Thus, Eqn. (2) can 
be modified into the form: 

x = W\<j — w 2 x — 4w 3 xy(l — kz), (3) 
y = Aw^a 2 — Aw 3 xy(l — ux), 
z = w 5 a — WqZ. 

We want Eqn. (3) as simple as possible. We include the reaction term 4w 3 nxyz in the 
first equation only because its effect in the second equation is assumed to be small. This 
assumption stems from a microscopic process: a pair of O a d s can not be blocked by a single 
inert molecule. The reason why we include 4w 3 ux 2 y in the second equation only is also 
taken from a microscopic fact: for the adsorption of CO, the presence of two free adjacent 
sites is not essential. 

We note that the two effects were also described by Jansen and Nieminen [10]; however, 
they did not include them explicitly in the rate equations. We also note that Eqn. (3) is 
much simpler than the rate equations in the pair approximation of Eqn. (1), which do not 
show oscillations in the parameter range of interest [11]. Eqn. (3) is an approximate model, 
and it is desireable to derive it systematically. 

IV. A HOPF BIFURCATION 

Because of its nonlinear terms, analytic treatment of Eqn. (3) is difficult. Nevertheless, 
it can be studied numerically. For example, one may use the software "XPPAUT" by B. 
Ermountrout [12]. 

We found that for the occurrence of a Hopf bifurcation it is not necessary that both k 
and u were nonzero in Eqn. (3). However, we found that parameter space for oscillations 
are larger if both k and u> are nonzero. 

A. case and 

In the present case, it is difficult to find analytic expressions for the steady states of Eqn. 
(3). In Fig. 3 we show a numerical integration of Eqn. (4). Here, in contrast to its dynamics 
shown in Fig. 1, the inert molecule concentration remains high. We suppose that the main 
reason for this is that w 3 is finite in Fig. 3, but w 3 = oo in Fig. I. Also, a better choice 
of the rate constants can be made in order to compare dynamics in the rate equations and 
Monte Carlo simulations [2]. 

A Hopf bifurcation in this system is shown in Fig. 4. Thus, the nonlinear effects we 
outlined in the previous section may indeed lead to stable oscillations. A two parameter 
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bifurcation diagram is shown in Fig. 5. We found that with the increase of W3, the parameter 
region for a Hopf bifurcation widens. We found that such an increase is accompanied by the 
emergence of saddle nodes. 



B. case and uj = 



In the present case, the steady states are given by, 



x = w 2 1 (w 1 z — 4:W 4 Z 2 + AKW4Z 3 ), 
y = l-x - 2z , 



(4) 



where zq is a solution to 



W2W4Z0 + 1^3 (u>l + Aw^Z^KZq — 1) 

{w 2 (2z - 1) + z (w 1 + 4w 4 z (kz - 1))) = 0. 



(5) 



Once again we were able to analyze stability of x , y and z numerically only. Oscillations 
for the present case are shown in Fig. 6. We note that in the present case, except the Hopf 
bifurcation, saddle node bifurcations occur in the system. 



As in the previous case, the steady states of Eqn. (3) can be expressed through a 
polynomial. Again, we were able to solve the polynomial numerically only. As a result, a 
linear stability analysis of these steady states has turned out to be difficult. Numerically 
detected oscillations for the present case are shown in Fig. 7. Note the exceptions of w^s 
value and w^O; these are the reaction rates used in Ref. [10] for Monte Carlo simulations. 
We note that for the parameters used in Fig. 6 - 7, a parameter domain for the Hopf 
bifurcation was smaller in the present case then it was in the previous case. 



A microscopic study of a reaction system based on Monte Carlo simulations can reveal 
the molecular mechanisms of catalytic processes. As an example, in this work, such study has 
allowed us to outline the essential nonlinear processes in site blocking caused oscillations. 
However, fora more complete study of a reaction system, the role of diffusion should be 
taken into account [13]. In that sense, Eqn. (4) can be a useful starting point for a 
reaction diffusion study of site blocking induced oscillations. The preliminary results of a 
reaction diffusion system in one dimension, with a diffusion of x, y, and z, show synchronous 
oscillations, stable long lived islands and diffusion induced chemical turbulence [14]. More 
detailed results on this, and its comparison to the Monte Carlo simulations, which include 
diffusion of adsorbents, will be reported elsewhere [15]. 

Finally, we suppose that the nonlinear effects and stable oscillations due to site blocking 
can be detected in experimental studies on CO oxidation. 



C. case k = and uj ^ 



V. DISCUSSIONS 
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FIGURES 



FIG. 1. Coverage oscillations on the square lattice. Open circles: 6co ada - Long dashed: 
0o ads , Solid lines: 0M ada - Rate constants: w\ = 0.44, u>2 = 0.005, ws = 00,11)4 = 0.56, 
w 5 = 0.004, w 6 = 0.004. Time is scaled to the total number of sites S = N 2 , N = 256. 

Q = £ Z) * = CO a d s , O a ds,M ads . 

FIG. 2. Snapshots of the surface at two different moments, a.) A phase without CO a d s 
coverage, b.) Emergence of an O a d s island on the CO a ds dominant surface. Parameters are the 
same as in Fig. 1. 

FIG. 3. Coverage oscillations in Eqn. (3). Open circles: x. Long dashed line: y. Solid line: z. 
Parameters are the same as in Fig. 1 except W3 = 5, u> = 1, k = 2. 

FIG. 4. A Hopf bifurcation diagram. Solid line represents a stable steady state. Circles repre- 
sent the Hopf bifurcation branch: filled circles show the supercritical Hopf bifurcation, open circles 
show unstable oscillations. Dashed and thin lines show an unstable steady state. Parameters are 
the same as in Fig. 3. 

FIG. 5. A two parameter bifurcation diagram. Oscillations occur in the area marked by osc. 
Other parameters are the same as in Fig. 3. 

FIG. 6. Oscillations in the case uj = 0. Open circles: x. Long dashed line: y. Solid line: z. 
Parameters are: w\ = 0.2, W2 = 0.01, W3 = 5, W4 = 1, W5 = w§ = 0.01 and k = 5.5 

FIG. 7. Oscillations in the case n = 0. Open circles: x. Long dashed line: y. Solid line: z. 
Parameters are: w\ = 1, W2 = 0.001, ^3 = 5, W4 = 0.52, = u)q = 0.0003, and to = 1. 
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Figure 2a-b 
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Figure 3 
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Figure 4 
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Figure 5 
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Figure 6 
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Figure 7 
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